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We discuss recent advances in tlie radiative-liydrodynamic modeling of core collapse super- 
novae in multi-dimensions. A number of earlier attempts at fully radiation-hydrodynamic mod- 
els utilized either the grey approximation to describe the neutrino distribution or utilized more 
sophisticated multigroup transport methods restricted to radial rays. In both cases these models 
have also neglected the 0(v/c) terms that couple the radiation and matter strongly in the op- 
tically thick regions of the collapsed core. In this paper we present some recent advances that 
resolve some shortcomings of earlier models. 



1. Introduction: The Supernova "Problem" 

For over a decade researchers have struggled to model the convective post-bounce epoch of core 
collapse supernovae. The radiative-hydrodynamic flows that occur in the region below the stalled 
prompt shock have held both promise and pitfall for the supernova modeler. The promise of this 
phenomenon is that it might explain the long sought-after mechanism that converts the core bounce 
into the observed explosion. The pitfalls are legion, mostly involving a complex convective flow 
structure that is three-dimensional in nature and couples neutrinos to matter strongly. For this reason 
there remain many open issues in modeling the convective epoch of core collapse supernovae. 

The supernova "problem" persists, despite more than four decades of concentrated research. 
The problem is this: We have no convincing explanation as to how the core collapse, which ends 
the evolution of a massive star, rebounds in such a way as to generate the explosion we observe in 
nature. The most realistic supernova models collapse and rebound, but create a shock wave that 
doesn't eject matter, either on the hydrodynamic timescale 10 ms), or the diffusive timescale of 
the escaping neutrino radiation 1-10 s), or on any other timescale we can model. 

This is not a problem of overall energetics. The gravitational energy released during core col- 
lapse and the subsequent neutron-star cooling phase is several factors of 10^^ erg. In contrast, the 
kinetic energy of the explosion required for consistency with observation is only 10^' erg. Instead of 
insufficient energy, the problem is one of energy conversion and transport — ^how a sufficient portion 
of the released gravitational energy is imparted to the material ejectus, giving it the requisite kinetic 
energy. 
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It has been understood for many years that neutrinos play a vital role in this process. In fact, 
essentially the entire remaining 99% of released energy (that which is not converted to kinetic energy 
of the matter) is radiated away as neutrinos. Thus, an accurate treatment of neutrino processes is a 
necessary component of any realistic model for a supernova. 

In this article, we present what we currently regard as the most important issues in core-collapse 
supernova modeling. In Sec. 2, we discuss the major components that need to be part of any serious 
modeling endeavor. In Sec. 3, we present an outline of V2D, our new two-dimensional (2-D) super- 
nova simulation code. Section 4 contains some preliminary results using this code. Our conclusions 
are in Sec. 5. 

2. The Components of a Supernova Simulation 

Broadly speaking, there are four main components to current supernova simulation models: (1) hy- 
drodynamics, to track the collapse, rebound, and ejection of stellar material, (2) neutrino transport, 
to track the production of neutrino radiation and to follow its propagation and emission from the 
star, (3) nuclear microphysics, to describe the diverse states of matter encountered throughout a sim- 
ulation, and (4) neutrino microphysics, to describe the reactions and interactions involving neutrinos 
and matter It must be stressed that all of these components are tightly coupled to one another Thus, 
the most effective models are designed with this couphng built in ab initio. 

Hydrodynamics. For simplicity of implementation, it has been customary for supernova codes 
to employ explicit hydrodynamics, with either a Newtonian or a general relativistic formulation. 
Implicit algorithms have usually been avoided since they require the computationally expensive 
solution of large systems of non-linear equations. 

Regardless of which approach is used, the hydrodynamic portion of the problem requires solu- 
tion of some form the following equations, expressed here in Newtonian formalism: 

^+V.(pv)=0 (1) 
— +V.(£v)+PV.v=-^y (3) 



dt 




V.(pw) + VP + pV4> + V.|^ I de{Pe + Ps) \=0. (4) 



Equation Q is the continuity equation for mass, where p is the mass density and v is the matter 
velocity, and where these quantities, and those in the following equations, are understood to be 
functions of position x and time t. Equation (|2} expresses the evolution of electric charge, where Ye 
is the ratio of the net number electrons over positrons to the total number of baryons. In the presence 
of weak interactions, the right hand side is non-zero to account for reactions where the number of 
electrons can change. Here, we express the net emissivity of a neutrino flavor (of energy e) and 
its antineutrino by §£ and Se, respectively. This expression is integrated over all neutrino energies 
and summed over all neutrino flavors /. The mean baryonic mass is given by m^. Evolution of the 



2 



internal energy of the matter is given by the gas-energy equation, Eq. (|3jl, where E is the matter 
internal energy density and P is the matter pressure. Again, the right hand side of this equation 
is non-zero whenever energy is transferred between matter and neutrino radiation as a result of 
weak interactions. We note that it is also possible to substitute for Eq. (|3j an expression for the 
evolution of the total matter energy (internal plus kinetic plus potential). Finally, Eq. (|4} expresses 
gas-momentum conservation, where <I> is the gravitational potential, and and Pg are radiation- 
pressure tensors for each energy and flavor of neutrino and its anti-neutrino, respectively. 

These equations must be discretized for solution within a computational framework. Tradition- 
ally, with one-dimensional models, it has been convenient to use Lagrangean methods, in which 
a computational mesh strictly co-moves with the mass elements of the fluid. With the advent of 
multi-dimensional models, however, it is common to use Eulerian hydrodynamics, where the mesh 
is fixed in an inertial frame of reference. This is because purely Lagrangean methods are difficult to 
implement in multi-dimensional schemes without the mesh suffering distortion and entanglement in 
convectively active regions. 

For all the benefits of Eulerian meshes, they also present a number of thorny issues. This is 
especially true for spherical polar meshes, the most natural choice for supernova modeling. The most 
obvious issue is the coordinate singularity that exists when the polar angle, 0. In addition, polar 
meshes exacerbate the problem of the timestep-restricting Courant-Friedrichs-Levy (CFL) condition 
at the center of the core. To deal with these issues, there are numerous resolutions and combinations 
of resolutions under active consideration. These include implicit methods, unstructured meshes, 
body-fitted meshes, and adaptive mesh refinement (AMR). 

As mentioned above, it is also necessary to choose between a total energy and an internal energy 
formulation . For the supernova problem, an internal energy formulation, as given in Eq. (|3j, is 
preferred. This is because much of the energy is internal, as opposed to kinetic. Solving the gas- 
energy equation helps insure an accurate calculation of the entropy, which is critical in degenerate 
regimes where a small change in energy can lead to a large change in temperature. 

The hydrodynamic algorithm must also have convergence properties that can deal with a realistic 
equation of state. This is particularly important in the regions of non-convex phase changes, such as 
the transition between nuclei and continuous nuclear matter. 

Neutrino Transport. This component is the most difficult to implement in a supernova model 
and the most time-consuming computationally. This is because supernova neutrinos cannot, in gen- 
eral, be described by an equilibrium distribution function. A solution requires a complete phase- 
space description of each neutrino's position and momentum. To obtain such a solution, one must 
solve the six-dimensional Boltzmann Transport Equation or some reasonable approximation thereof. 
This extra dimensionality easily leads to the transport calculation completely dominating a simula- 
tion in terms of computer memory, execution time, and I/O requirements. 

The Boltzmann Transport Equation (BTE) can be expressed in terms of the radiation intensity, 
/ — /(e,x,£2,f ), where e is the energy of a neutrino, x its position, and £1 the solid angle into which 
the neutrino radiation is directed. In terms of /, the Newtonian BTE can be expressed as 



where a, is the / ' component of the matter acceleration and pi the / ' component of the momentum 
of the neutrino. The right hand side of Eq. (|5} lumps together the contributions from all interactions 
that a neutrino might experience and is collectively referred to as the collision integral. 

A storm of issues faces one who implements a neutrino transport algorithm. Mezzacappa and 




(5) 
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BruenJiSlhave the only "full" solution to the BTE implemented in supernova simulations and then 
only with one-dimensional hydrodynamics. Upon moving to multi-dimensional models, the full 
solution of the BTE becomes yet more challenging to implement and more time-consuming to com- 
pute. However, this is the way that the field must ultimately go. (Livne and colleagues purpor^lto 
implement a two-dimensional S,, solution of the BTE. However, this solution omits critical matter- 
radiation coupling terms and no numerical details of the method have been disclosed.) 

In the meantime, a number of approximate transport moment methods have emerged. The most 
successful of these is use of a finite series of angular moments of the BTE. When this approach is 

taken, a limiting scheme is then required to close the resulting equations. Breunn^ implemented 

7 1 

a Pi scheme, with flux-limiting. Myra et al. closed the zeroth angular moment^ of the BTE 



by implementing the Levermore and Pomraning flux limitepSI Bowers and Wilsor^J also used a 
flux-limiting scheme of their own device. More recently, Rampp and Janka^ have implemented 
a variable-Eddington-factor approach to solve the first two angular moments of the BTE. In all the 
above cases, however, the implementation has been made in only one spatial dimension. (Janka et 
alU^ are developing a two-dimensional Boltzmann transport code (MuDBaTH), but the numerical 
details are as yet unpublished.) 

Since all the schemes noted so far derive monochromatic transport equations, yielding a separate 
equation for each neutrino energy, they are referred to as multi-group schemes. Those that combine 
multi-group and flux-limiting are known as multi-group flux-limited diffusion (MGFLD) schemes. 

A much simpler alternative to multi-group schemes is so-call "grey" transport, which is derived 
by integrating the BTE over both neutrino energy and angle. To perform these integrals, one must 
assume a spectral shape for the neutrino distribution. Typically, this requires defining an arbitrary 
neutrino "temperature," and assuming that neutrinos can be parameterized by some kind of equi- 
librium Fermi-Dirac distribution. Among relatively recent models, this was first imp lemented by 
Cooperstein, van den Horn, and Baron,dand later by Swest}l22land by Herant et alU^This approx- 
imation is still in active use by the latter group. ^ 

Grey schemes have numerous shortcomings. First, work with multi-group schemes has shown 
that in areas where accurate neutrino transport is critical, neutrinos do not assume any kind of distri- 
bution that can parameterized once and for all as required by grey transport. Spectral distributions 
constantly evolve and, thus, a multi-group description is required to obtain even a qualitatively cor- 
rect description. More troubling, Swesty^ has shown that by adjusting the grey parameterization 
within very small bounds, it is possible to "dial" an explosion (or failed explosion) with the appro- 
priate choices of these unknowable and unphysical tuning parameters. Hence, although grey codes 
have utility for making a sweeping exploration of parameter space, any scientific conclusions that 
rely on them should be viewed as highly suspicious, and not regarded as in any way definitive. 

Regardless of which transport scheme is implemented, another critical issue that must be faced 
is matter-radiation coupling. Coupling occurs in Eqs. (|2}-0 for lepton number, energy, and mo- 
mentum evolution. Coupling that occurs on the right-hand side of these equations has a conceptu- 
ally simple analytic structure. However, momentum transfer is more troublesome in approximate 
schemes. This is because the radiation momentum equation is often truncated in a way that makes 
the accuracy of the calculated momentum transfer less certain. 

Matter-radiation coupling also enters implicitly in the neutrino transport equation through spec- 
tral rearrangement terms and in the dynamic diffusion term. Both these terms are frequently and 
erroneously neglected in supernova models, even though the dynamic diffusion term is the leading 
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order contribution in optically-thick regions. 

Equation of State. Adequate modeling of stellar-core collapse requires an equation of state 
(EOS) that handles a density range of roughly 10^-10^^ g cm^^, a temperature range of 0. 1-25 MeV, 
and an electron-fraction range of 0.0-0.5. The EOS must also be able to handle different regimes 
of equilibrium states. Throughout most of the core, the material is in nuclear statistical equilibrium 
(NSE) and is usually modeled by one of the NSE equations of state. Although the gross features of 
nuclear matter are thought to be well-understood, there is still much open ground for investigation. 
Fertile regimes for such work include the EOS at supernuclear densities. In addition, little is known 
about the nature of nuclei at subnuclear densities when Ye is small. 

Matter in the silcon shell and beyond does not attain NSE until the bounce-shock wave passes 
through it. Dealing with the transition between NSE and non-NSE EOS's and with the network of 
nuclear reactions that joins them is a challenge that is only beginning to be addressed.-H- 

Neutrino Microphysics. Since the energetics of a core-collapse supernova is primarily a neu- 
trino phenomenon, it is necessary to have correct opacities and rates for the various neutrino pro- 
cesses that are important. The collection of reactions that are important, or possibly important, to 
the supernova problem is rich and has evolved through the years. Arguably the most important de- 
velopment came with the discovery of weak neutral currents, from which it could be inferred that 
the dominant contribution to neutrino opacity in a collapsing stellar core is from coherent elastic 
scattering of neutrinos from nuclei.El 

The list of possible neutrino interactions is nearly endless, but those of demonstrated importance 
include the coherent scattering just mentioned, as well as conservative scattering from free nucle- 
ons. Also of undisputed importance are electron capture by protons (and protons bound in nuclei), 
neutrino production through electron-positron pair annihilation, and neutrino-electron scattering. 

In recent years, with the experimental evidence pointing strongly to the existence of neutrino 
oscillations, it is also important to investigate the possible role of flavor-changing interactions to the 
supernova problem. Investigation into this has begun,E^but has yet to be incorporated in a detailed 
simulation. 

3. V2D: A New Code for Two-Dimensional Radiation Hydrodynamics 

Our new radiation-hydrodynamic simulation code, V2D, is a two-dimensional, Newtonian, pure Eu- 
lerian, st aggered-m esh code based on a modified version of the algorithm for ZEUS-2D by Stone and 
Norman-ESEm^Following Stone and Norman, it has been designed for use in a general orthogonal 
two-dimensional geometry, which makes its utility extend beyond the supernova problem. 

V2D is an entirely new implementation, coded according to the Fortran 95 standard. It is a 
distributed-memory parallel code that uses calls to MPI-1 for message passing between processes. 
It has been designed for easy portability between computing platforms and currently runs on systems 
ranging from as small as a Linux-based laptop to as many as 2048 processors of an IBM SP. To aid in 
this portability, the input and output is formatted using parallel HDF5, which is built on the MPI-I/O 
portion of the MPI-2 standard. 

One of the major design goals of V2D is componentization and, to adhere to this principle, we 
insist on completely separating microphysics from the numerical implementation of our radiation- 
hydrodynamics algorithm. This isolation of mathematics and computational science from physics 
has allowed significant contributions from applied mathematicians to enhancing the performance of 
our code. 
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The V2D algorithm rehes on operator sphtting, with advection steps spHt from source-term 
steps. Hydrodynamic and neutrino-transport source-term steps and coupHng are interleaved. At 
the start of each simulation timestep, the gravitational mass interior to each point is calculated. 
Since the collapsed core is nearly spherically symmetric and highly condensed, we approximate the 
gravitational mass assuming that mass interior to the point of interest is in a spherically symmetric 
distribution. In future versions of our code, we will implement a more accurate Poisson solver to 
calculate the (slightly) non-spherical gravitational potential. 

In a break from Stone and Norman's method, V2D next performs the advection sweeps in the 
radiation-hydrodynamic quantities (mass density, matter internal energy, velocities and momenta, 
electron fraction, and neutrino distributions). Following this, a neutrino transport step is performed 
for each flavor and each matter-radiation energy exchange is calculated. 

The matter pressure is next updated, upon which the gravitational and matter- and radiation- 
pressure forces are applied to the matter. Artificial viscosity is calculated next and its contributions 
applied to the fluid. Finally, the gas-energy equation is solved. 

This procedure is repeated for each timestep in a simulation, with the provision that advection 
sweeps are ordered alternately according to timestep (i.e., xi-direction first, followed by X2, or vice 
versa). 



3.1. Neutrino Transport Implementation 

As an extension of earlier work by usj ^^ l '^^ lwe implement neutrino transport by taking the zeroth 
angular moment of the BTE to yield the following neutrino monochromatic energy equation in the 
co-moving frame: 

^+V.(£,v)+V.F,-e^(P,:Vv)=§„ (6) 

where is the neutrino energy density per unit energy interval at position x and time f , Fg is the 
neutrino energy flux per unit energy interval, and and Sg are as defined earlier. The expression 
Pg : Vv indicates contraction in both indices of the second-rank tensors Pg and Vv. There is a 
corresponding equation to describe the antineutrinos. This pair of equations is repeated for each 
neutrino energy e, and neutrino flavor. We currently track electronic, muonic, and tauonic neutrinos. 

Equation (|6jl is closed using Levermore and Pomraning's prescription for flux-limited 
diffusion,II3 which allows us to express Fg as 

Fg = -DgV£g, (7) 

where Dg is a "variable" diffusion coefficient that varies in such a way as to yield the correct fluxes 
for the diffusion and free streaming limits and an approximate solution in the intermediate regime. 
This prescription also provides the elements of the radiation-pressure tensor Pg. 

Presently, we emplo y the same prescriptions for electron capture and conservative scattering that 
we have used in the pasti ^^ l ^ ^The rates for these process are calculated on the fly within the course of 
a simulation. Our model also implements neutrino production via electron-positron pair annihilation, 
as in Yueh and Buchlei-^ and Bruenn . Neutrino-electron scattering has been implemented, but is 
not currently turned on in the preliminary results we present here. These latter two sets of processes 
use tables of precomputed rates, which are interpolated via a tri-linear interpolation scheme over 
neutrino energy £, temperature T, and electron chemical potential jie. 

Since the neutrino CFL restriction on a transport timestep is far too restrictive to permit an 
explicit solution, we use a purely implicit method to solve the transport. The equations compris- 
ing the description of each neutrino-antineutrino species are assembled in matrix form. We note 
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that the second (advective) term in Eq. (|6} is omitted from this process since it has been already 
treated during the operator sphtting of the advective step described above. Blocking terms arising 
from Fermi-Dirac statistical restrictions on final neutrino states make this a system of non-linear 
equations. Fortunately, the system is sparse, which makes it amenable to solution by sparse iterative 
methods. A nested procedure is used, employing Newton-Krylov methods^. In the innermost loop, 
a linearized system is solved using preconditioned Krylov-subspace methods. The outer loop uses 
a Newton-Raphson iterative scheme to resolve the non-linearity of the system. Besides being an 
effective general procedure for sparse systems, our implementation of parallel preconditioners also 
insures that is amenable to large-scale solution on parallel architectures. This is the chief reason our 
code exhibits its high degree of scalability across many platforms. 

3.2. Equation of State 

V2D is designed to use an arbitrary equation of state and we use several in the course of testing the 
code. For production runs, however, we use the Lattimer-Swesty EOS^^ in tabular form. The 
thermodynamic quantities are tabulated in terms of independent variables, density, p, temperature, 
T, and electron fraction. Ye- We have tabulated this EOS in a thermodynamically consistent way 
according to the prescription in Swest> (We refer to this combination collectively as LS-TCT.) 
We note that although the Lattimer-Swesty EOS is commonly used, and tabulations of it are also 
common, most tabulations are not constructed in such a way as to guarantee thermodynamic consis- 
tency. When non-thermodynamically-consistent tables are used, spurious entropy can be generated 
or lost. Such problems have been sometimes incorrectly attributed to the Lattimer-Swesty EOS, 
rather than erroneous tabulation of the otherwise consistent EOS. 

To guarantee tabular thermodynamic consistency, LS-TCT uses bi-quintic Hermite interpolation 
in the Helmholtz free energy F , as a function of T and p. Functional dependence on Ye changes 
slowly enough to permit linear interpolation. This procedure is required since satisfaction of the 
Maxwell relations requires consistency among the second derivatives of F . In addition, we desire 
fidelity of the interpolation and continuity of derivatives to the underlying tabular data for both F 
and its derivatives, {dF /dT)p and [dF /dp)T- 

Apart from thermodynamic considerations, we also want to insure that there are no discontinu- 
ities that might cause difficulties in the hydrodynamics. Hence, LS-TCT also insists that interpola- 
tions of the second derivatives of F approach the correct tabulated values. (This is the equivalent 
of saying that we require the derivatives of pressure and internal energy with respect to temperature 
and density be continuous.) 

A final requirement is that we wish the interpolation function and its first and second derivatives 
be continuous across table-cell boundaries. This insures that nothing untoward happens as a fluid 
element migrates from one thermodynamic regime of interest to another. 

3.3. Software Validation and Verification 

When engaged in a major project, such as V2D, it is important for software developers to be con- 
stantly vigilant about the quality of software being produced. It is important that the software meet 
the requirements that it is intended to address (validation) and that the code yields correct answers 
(verification). 

We take these issues seriously and, in an effort to address them, have implemented strict source- 
code control and testing procedures to ensure that our software meets our rigorous standards. One 
of the most important elements of our program has been the implementation of a suite of regression 
tests, which we currently run four times daily. This is not a static suite, but is constantly growing. 
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Our eventual aim is to cover every major element of code in V2D. Currently our suite consists 
of about two dozen separate problems that include tests of the hydrodynamics, neutrino transport, 
parallel solvers, message passing, and parallel I/O. Wherever possible, we try to include problems 
with analytic or at least verifiable solutions. 

Although implementation of these procedures is labor intensive and time-consuming, we feel it 
is a justifiable investment. With current regression tests, our procedures have already been effective 
in finding errors in our code. They have also served as an effective safeguard against introducing 
new errors as we continually enhance V2D's functionality. 

4. Initial Results 

We are using V2D to carry out our first 2-D multigroup models of the post-bounce epoch. To date, 
simulations have not reached a sufficient time that would allow us to be certain about whether an 
explosion is obtained. Nevertheless, the results warrant some discussion as they reveal important 
features of the post-bounce epoch. 

4.1. Initial Model 

For a progenitor model we employ the widely-used Woosley and Weavei^^l S 15S7B2 ISM© pro- 
genitor Much previous work has focused on the evolution of this progenitor through collapse, core 
bounce, and convective phases. The Fe core, the Si shell, and a portion of the O shell area are zoned 
into a 256 radial-mass-zone mesh with zoning that is tuned (by trial and error) so as to yield a high 
spatial resolution grid in the proto-neutron star and the inner 200 km of the collapsed core at bounce. 
This tuned zoning sets up a radial grid that is compatible with subsequent 2-D Eulerian simulations. 
The neutrino-energy spectrum, ranging from 0-375 MeV, is discretized into 20 energy groups with 
group widths that increase geometrically with energy so as to resolve accurately the Fermi surface 
of the electrons and neutrinos in the proto-neutron star. The initial values for T, p, and Ye, are 
interpolated from the original S15S7B2 data onto the Langragean mass grid and the initial radial co- 
ordinates of each mass shell are computed consistently with density. The neutrino energy densities 
Ev are initiaUzed to a small non-zero value that yield an initial neutrino luminosity that is many or- 
ders of magnitude below what the precollapse thermal pair-production luminosity would be. When 
the simulation is started in its pre-collapse quasi-static phase, the luminosity stabilizes within a light 
crossing time (5-10 ms). 

4.2. Lagrangean Collapse Calculations 

The initial model is collapsed using a 1-D Newtonian Lagrangean radiation hydrodynamics code 
RHID that uses the mass and energy meshes described above. The evolution algorithm for each 
timestep utilizes operator splitting to first carry out a Lagrangean hydrodynamics step followed by 
Lagrangean neutrino evolution steps for each of the three neutrino flavors. After each Lgrangean 
neutrino evolution step the matter internal energy and electron fraction are corrected for any energy 
and lepton number exchange that has occurred. 

The model contains neutrino microphysics as described in BruenrPI with two exceptions. The 
neutrino-nuclei scattering opacity has been modified to take into account the form-factor introduced 
by Burrows, Mazurek, and Lattimei^ and we have neglected the effect of neutrino/anti-neutrino 
annihilation. Full neutrino-electron scattering as described in this paper is included in the code but 
is not turned on in the model described in this paper. Additional effects such as nucleon recoil, 
ion-ion correlations, etc. are being considered as follow-ons to this baseline model. 
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The EOS utilized is the TCT tabularized version of the nuclear EOS Lattimer-Swesty (LS-TCT) 
with the = 180 MeV parameter set, with the exception of the electron EOS that has been updated 
to improve accuracy and the range of applicability. The original LS EOS chose a value for the alpha 
particle binding energy Ba =28.3 MeV that did not correctly account for the neutron-proton mass 
difference. This parameter has been subsequently corrected but in the model presented in this paper, 
the original value has been retained so as to allow comparison to other work. 

Using RHID with the microphysics and meshes described above the core collapses in ap- 
proximately 233 ms. The central conditions at bounce are approximately T « 10.3 MeV, p « 
2.71 X lO'"* g cm^^. Ye « 0.306, and « 0.37. Thes e condition s are in good agreement with those 
obtained in Lagrangean MGFLD and MGBT models EEIEZini. 

4.3. Eulerian 2-D Calculations 

Our 2-D models are carried out in spherical polar coordinates. The initial conditions for our 2- 
D simulations are taken from the 1-D Langrangean simulations as the central density of the core 
reaches the nuclear saturation density. The T, p, and Ye profiles at this point are shown in Figure 
[2 and are taken at approximate time of 233.087 ms. The choice of this epoch to begin our 2- 
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Figure 1 . Initial radial profile for 2-D simulation. Solid line indicates temperature in units of MeV. Dashed line indicates 
iO xYg. Dotted line indicates logiQ(p) in units of g cm^'. 

D models was made for two reasons. First, the prompt shock propagates and stalls into a quasi- 
static equilibrium on the Eulerian mesh. We have found that if we attempt the transition from 
a Lagrangean algorithm to an Eulerian algorithm at the point where the shock has stalled, there 
will be noticeable differences in the quasi-static equiUbrium caused by the presence of the nuclear 
"pasta" phase transition. These differences in the equilibrium point can set off spurious unphysical 
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shocks. By choosing to let the prompt shock propagate a stall on an Eulerian mesh, we avoid this 
problem. The second reason for choosing the initial point for our 2-D models near the moment of 
bounce is that the radial coordinates of the zones in the outer part have achieved desirable values for 
an Eulerian simulation. 

The radial zoning for the 2-D models is taken directly from the radial coordinates of the 1-D 
Langrangean model zones. In this way we avoid any need for remapping of data between radial 
zones. This eliminates the introduction of any spurious forces into the initial conditions for the 
2-D models. The initial data for T, p, Yg, v^, and the neutrino energy densities Ey and By for 
each neutrino flavor are taken directly from the 1-D Langrangean model. The initial velocity in 
the 9 direction is set to zero. The data are mapped in the polar-angular direction in a spherically 
symmetric fashion. The angular grid consists of 256 zones uniformly spaced in angle over the range 
of < < ;r. We place a small sinusoidal perturbation in the electron fraction Ye to seed convection 
in the region between 100 and 200 km of the form (}'e)pertuib — (ie)Lagiangean +CpSin(40) where 
Cp = 10^^. The energy group structure is also left unchanged from the 1-D Langrangean runs so 
there is no need to remap the data in the energy dimension. In the remainder of this paper we shall 
refer to this model as Production Run 37 (PR37). 

The use of spherical polar coordinates in combination with explicit hydrodynamic algorithms to 
model spatial domains that include the origin gives rise to a numerical stability problem. At r = the 
spherical coordinate system is degenerate and all zones that include the origin as a vertex should be in 
instantaneous sonic communication with one another However, standard explicit numerical finite- 
difference, finite-volume, or finite-element techniques are limited to nearest-neighbor type spatial 
coupling and do not include numerical coupling between all zones containing a given vertex. In order 
to circumvent this problem we numerically introduce "baffles" into the center of the collapsed core, 
as though it were a tank of fluid, to prevent movement of fluid in the angular direction inside a certain 
radius. Since no fluid movement occurs in the 9 direction inside the baffle radius, there is no CFL 
restriction based on the zone size in the 9 direction for zones inside that radius. Nevertheless, the 
zones on either side of the baffle are sonically connected as sound waves flow around the outer edge 
of each baffle. We strive to keep the baffle radius small, so that the flow in the 9 direction remains 
unimpeded in any region where convective instabilities may develop. For the mode described in 
this paper the baffle radius is approximately 8.5 km which yields an average CFL timestep of about 
5 X 10^^ s. As we will see, this baffle radius is well inside of any proto-neutron star (PNS) instability 
region. 

The 2-D calculations have been carried out on the IBM-SP system at the National Energy Re- 
search Scientific Computing Center (NERSC). The models are run on 1024 processors with paral- 
lelism handled via message passing via calls to MPI libraries. Model PR37 required approximately 
50,000 processor-hours of CPU time to reach a simulation time of 16 ms post-bounce. 



4.4. The Onset of Convection 

The 2-D models are carried out from the point of bounce. As expected the prompt shock weakens 
while propagating outward and finally stalls. In the 2-D models, the radius of the shock at 5 ms 
post-bounce is near 60 km and propagating outwards very slowly. This is in good agreement with 
our 1-D Lagrangean models. 

One difference that we see from previous works including the grey models carried out by 
SwestjISH is that convective instabilities are born much earlier when the 2-D models are initial- 
ized near the point of core bounce. By 10 ms after bounce, the model has developed two separate 
unstable layers. The outer layer, shown in Fig. (|2j, is the classic entropy driven Rayleigh-Taylor 
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convection. 



The inner layer, shown in Fig. (|3}, seems to be an instabihty in the outer lay ers of the proto- 
neutron star similar to those seen in 2-D simulations carried out by by Keil et a/.li^land Mezzacappa 
etalW. 



entropy Time = 242.411 msec 




-60 -40 -20 20 40 60 80 

Kilometers 



Figure 2. The entropy per baryon for model PR37 at approximately 9 ms after core bounce. 



T here is c ontroversy about the existence of PNS instabilities. Originally, the work by Wilson and 
Maylj22III] claimed to find doubly-dif fusive instabilities in the region below the neutrinosphere. 

Later work by Breunn and collaborators'^'"^'^' cast doubt on the existence of such phenomena. The 

1 '\ 

simulations of Keil, which utilized the grey flux-limited diffusion approximation to transport neu- 
trinos along radial rays, found a PNS instability that grew over the time period of approximately one 
second to encompass the entire proto-neutron star In contrast, the subsequent work of Mezzacappa 
et al. found a PNS instability that quickly damped out within a short time. The Mezzacappa et al. 
simulation utilized a multigroup flux-limited Pi approximation to transport neutrinos along radial 
rays outward from an inner radius of r=20 km. The inner boundary conditions in this simulation 
were established as time-dependent data from the 1-D Langrangean code of BruenrH It is impor- 
tant to note that neither the simulations of Keil et al)-^ or Mezzacappa et a/.l^ took into account 
the fully radiation-hydrodynamic coupling via the compression and dynamic diffusion terms that 
are present in Eq. (|6j. Our simulations have confirmed the expected result that the effects of the 
dynamic diffusion term dominate the radiative diffusion term in the regions in which the optical 
depth is large. 

Figure |3l shows the instantaneous structure of the velocity field in model PR37 at the same time 
as the data shown in Fig. (|2j- The velocity field is illustrated by means of a Lagrangean-Eulerian 
Advection (LEA) visualization technique that shows the direction of the vectors as streaks. One can 
clearly see eddies associated with the PNS instability layer at a radius of about 20-25 km. This is 
well outside the baffle radius of 8.5 km. In fact, there are approximately a minimum of 40 radial 
zones separating the innermost of the vortices and the outer edge of the baffles. We do not believe 
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flW 3JJ 6.Sf lO.Q 133 >S.7 10.0 333 5*7 



Figure 3. The velocity structure of tlie PNS instability layer in model PR37 at t=242.411 ms. The image depicts the 
instantaneous structure of the velocity tield by means of a texture map visualization technique known as Lagrangean-Eulerian 
Advection (LEA). 



that the baffles in any way impede the dynamics of the vortices. Nevertheless other simulations are 
underway where the baffle radius is made substantially smaller to verify this claim. We are also 
developing an implicit hydrodynamic algorithm that will avoid the need for baffles altogether. 

Whether the PNS instability will grow or diminish with time is as yet unclear since we have 
only evolved model PR37 to a time of approximately 16 ms at the time of this writing. The velocity 
structure at that time is shown in Fig. (0}, which clearly reveals that much coherency has been 
lost in the vorticial structure of the PNS instability layer This seems indicative of the decay of 
this PNS instability in a fashion similar to that described by Mezzacappa et al. . However, it is 
necessary to evolve this simulation substantially farther in time before any definitive statements can 
be made about the long-term behavior of this sector of the proto-neutron star. During this relatively 
short timescale the outer convective zone seen in Fig. Q does not exhibit significant growth. With 
evolution to the next 30 ms, we should be able to make comparative statements regarding model 
PR37 and earlier grey models carried out by SwestjI^H 



5. Conclusions and Future Directions 

The issues and work described in this paper fall far short of offering complete coverage of the active 
issues that remain in the area of the explosion mechanism of core collapse supernovae. Indeed, we 
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have ignored many important issues such as magnetic fields, rotation, and neutrino flavor mixing. 
Clearly there is a large gulf of unexplored physics incorporated in those subjects. 

Our own future efforts will be focused in two areas in the near future. The first of these is 
understanding the effects of numerous microphysics enchancements that will be added to the models. 
The second of these efforts involves extending the models to 3-D where a more realistic convective 
flow structure can arise. 
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